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We study the effects of demographic stochasticity on the long-term dynamics of biological co- 
evolution models of community assembly. The noise is induced in order to check the validity of 
deterministic population dynamics. While mutualistic communities show little dependence on the 
stochastic population fluctuations, predator-prey models show strong dependence on the stochastic- 
ity, indicating the relevance of the finiteness of the populations. For a predator-prey model, the noise 
causes drastic decreases in diversity and total population size. The communities that emerge under 
influence of the noise consist of species strongly coupled with each other and have stronger linear 
stability around the fixed-point populations than the corresponding noiseless model. The dynamics 
on evolutionary time scales for the predator-prey model are also altered by the noise. Approximate 
1/f fluctuations are observed with noise, while 1/f^ fluctuations are found for the model without 
demographic noise. 
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I. INTRODUCTION 



Noise may be a relevant perturbation to many kinds 
of population dynamics. Effects of population fluctua- 
tions have been investigated for several models, for exam- 
ple, predator-prey models P, [2| , epidemic models M , the 
Ricker model [j] , evolutionary game theories j5l-[lG| , and 
pattern formation [TT| - fl3j . Since the birth-death process 
of individuals is stochastic, the population of each species 
always fluctuates due to the flniteness of the number of 
individuals, known as demographic stochasticity. Demo- 
graphic stochasticity is an endogenous phenomenon, and 
the species populations may fluctuate, even in a constant 
environment. Population dynamics with demographic 
stochasticity are more realistic than the corresponding 
deterministic description which is valid for infinite pop- 
ulations in constant environments, and they often show 
nontrivial dynamics which cannot be predicted by the de- 
terministic equations. For a particular predator-prey sys- 
tem the demographic stochasticity causes oscillations, 
while the system is asymptotically stable under the corre- 
sponding deterministic population dynamics. Since the 
demographic stochasticity effectively adds uncorrelated 
noise to the population dynamics, the oscillations at the 
eigenfrequencies are amplified by a large factor. Hence, 
the effect of the demographic stochasticity is much larger 
than the one estimated by naive 0(1/ Vn) estimates. For 
a neutrally stable system @, the noise effect becomes 
even more drastic: only one species can survive and the 
others die out after a sufficiently long time, while the cor- 
responding deterministic model predicts the coexistence 
of the species, with regular oscillations. Stochasticity 
may also influence the outcome of the evolutionary dy- 



namics. In small populations, the evolutionary branch- 
ing is delayed compared to the case of larger populations, 
and the delay strongly depends on the absolute popula- 
tion size Several empirical data sets are also com- 
pared with theoretical models and are described better 
by models with stochastic population dynamics [H, 
Thus, population fluctuations, which inevitably exist in 
any finite system, may drastically alter the predictions of 
deterministic models and often cause decreases in biodi- 
versity. A major goal of this paper is to investigate the 
effects of demographic stochasticity in models of biolog- 
ical community assembly on evolutionary time scales. 

Several models to bridge ecological and evolutionary 
time scales have been suggested, such as the tangled- 
nature model ^M^, simplified versions of that model 



18H22i, the Webworld model 23H26i, the scale-invariant 



model [27[ , and others [28|,|29| . More concretely speaking, 
these are population dynamics models with additional 
rules for the introduction and extinction of species. New 
species, whose interaction coefficients are assigned by a 
rule, are added to the community at a certain rate; and 
the extinction of resident species can happen due to the 
population dynamics. Evolution is modeled by repeating 
the introductions and the extinctions of species in these 
models. Potential numbers of species are much larger 
than the number of species coexisting at the same time. 
If the population dynamics have a noise term, the emer- 
gent communities can be nontrivially different from the 
communities selected without this noise. This is because 
the speciation events that trigger large changes at the 
population level invariably involve single or very small 
numbers of individuals that are highly susceptible to sta- 
tistical ffuctuations 14 1 . The main issues we address here 
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arc the noise effects on (i) the properties of the emergent 
communities and (ii) the statistics of the evolution dy- 
namics, especially the intermittency during evolution. 

In this paper, we use the simplified versions of the 
tangled-nature models and study the effects of demo- 
graphic stochasticity. For these models, the fixed point 
and the linear stability around it for a given commu- 
nity are analytically obtainable. This helps us esti- 
mate the relevance of the noise in the population dy- 
namics. Furthermore, long-term dynamics on evolution- 
ary time scales are studied extensively, and 1// fluctu- 
ations and power-law duration distributions are found 
for these models. We show that the dynamics observed 
for the individual-based models may undergo qualitative 
changes from the corresponding models with determinis- 
tic population dynamics. 

The organization of this paper is as follows. In sec- 
tion ini the models are defined, and some topics related 
to these models are discussed. Results of the simula- 
tions are shown in section Hill In section HlI Al we show 
how population fiuctuations affect the diversity, and in 
section IIIIBl we explore the dynamics on evolutionary 
time scales. Section IIVI is devoted to a summary and 
discussion. Some mathematical details are discussed in 
Appendices EllEl 



II. MODELS 

The models considered here are extensions of the 
tangled-nature model studied in [l8l^,20] . The tangled- 
nature model is an individual-based model, originally in- 
troduced by Hall and co-workers and later simplified 
by Rikvold and Zia [1^. In the simplified models (Tsl - 
[23 . [30| , the population evolves stochastically in discrete, 
non-overlapping generations. In these models, each in- 
dividual of species / gives rise to F offspring with a re- 
production probability P/ before it dies. Otherwise it 
dies without offspring. The reproduction probability Pj 
for an individual of species I in generation t depends on 
the individual's ability to utilize the amount R of avail- 
able external resources, and on its interactions with the 
population sizes nj{t) of all the species present in the 
community at that time. The form of Pi is discussed in 
the following subsection. 

In the individual-based models, species populations 
evolve stochastically. The probability p/(/c) that k out of 
n individuals of species / succeed in producing offspring 
is given by the binomial distribution, 

Pdk)^(^fjPni~Pi)"-'', (1) 

where (^) is the binomial coefficient and Pj is the proba- 
bility that an individual of species I gives rise to offspring 
in that generation. Thus, the mean and variance of the 
number of offspring are nP/ and nP/(l — P/), respec- 
tively. In this paper, we consider an approximation to 



the binomial distribution by the Gaussian distribution 
with mean nPj and variance nP/(l — Pi) in order to 
control the strength of the stochasticity. The following 
stochastic difference equation is used for the population 
updates: 

nj{t + 1) = F[Pjnj{t) + K^ni{t)Pi{l- Pi)i{t)] , (2) 

where k and ^{t) are a control parameter for the noise 
strength and a Gaussian random number with mean 
and variance 1, respectively. When k = 1, this update 
algorithm is a good approximation for the corresponding 
individual-based model, while it is deterministic when 
K = 0. Although there is no easy interpretation except 
for K = and 1, we use several intermediate values of 
K in order to investigate the crossover between deter- 
ministic and individual-based models. The population 
size ni(t) is a positive real number while it is a posi- 
tive integer in the original individual-based models. The 
approximation by the Gaussian distribution to the bi- 
nomial one is known to be quite good when nj is suf- 
ficiently large. Typically when nPj and n(l — Pj) are 
greater than five, the approximation is good. Even when 
nj is small, we expect that the approximate model still 
captures the essence of the population fluctuations in the 
individual-based model, although the population dynam- 
ics for species with very small populations are critical for 
the emergence or extinction of species. Since it is not 
necessary to draw random numbers for every individual, 
this update rule is computationally more efficient than 
the true individual-based model and enables us to run 
simulations for longer times. 

It is straightforward to extend the model so that the 
number F of offspring per individual follows a stochastic 
process. (See Appendix 1X1 ) In that case, the fluctuations 
are even more enhanced and the difference from the de- 
terministic models are more important. In this paper, 
however, we limit ourselves to the case that F is fixed 
for simplicity. Even with this model, the differences be- 
tween stochastic and deterministic population dynamics 
arc observed as shown later. 

To mimic an evolutionary process, speciation and ex- 
tinction of species are introduced. New species are added 
to the system by "mutation" of resident species. These 
rules are formulated in section III Bl Extinction of a 
species happens when its population becomes less than a 
threshold value, nthr = 0.5. When the Jth species goes 
extinct, this species is eliminated from the system, i.e., 
the number of degrees of freedom decreases by one. The 
community configuration reorganizes as a result of the 
appearance and extinction of species. 

A. Reproduction probability 

As in the original models [2^, the reproduction prob- 
ability Pj is taken as 
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where 



Ai{R,{njit)}) = -bj + 



Ntotit) ^ A^totW No ■ 

(4) 

Here bj is the cost of reproduction for species / (always 
positive), and 77/ is the abihty of individuals of species / 
to utilize the external resource R. The interaction ma- 
trix M defines the interactions between species. The to- 
tal population size is denoted by Ntot{t) = 'Y^jnj{t), 
and A'o is an environmental carrying capacity that pre- 
vents A'tot (0 from diverging to infinity. The reproduction 
probability is a monotonically increasing 

function of A/, ranging over (0,1). Thus A/ is a measure 
of the fitness of species /. For a large positive A/ (small 
birth cost, strong coupling to the external resource, and 
more prey than predators). Pi goes to one and the popu- 
lation of species / increases. In the opposite limit of large 
negative A/ (large birth cost, weak or no coupling to the 
external resources, and/or more predators than prey). Pi 
goes to zero and the population decreases rapidly. The 
nonlinear dependence of Pi on A/ thus limits the growth 
rate of the population size, even under extremely favor- 
able conditions for species /. 

Two types of reproduction probabilities arc considered 
in this paper: Model A and Model B. Model A has no re- 
striction on the form of the interaction matrix M. There- 
fore, each species makes various types of interactions with 
others, including predator-prey, mutualistic, and compet- 
itive interactions. In contrast, the interspecies interac- 
tions arc limited to predator-prey interactions in Model 
B. This is realized by the limitation that the off-diagonal 
part of M must be antisymmetric {Mij ~ —Mji). Thus, 
if Af/j > and Mji < 0, then species / is the predator 
(or parasite) and J the prey (or host), and vice versa. 
Model A has a more general form, while Model B focuses 
on the energy transport via the foodweb. 

Model A was introduced and studied in [lO,!!^]. In this 
model, the reproduction cost &/ and the external resource 
R are zero; thus the first and the second terms of Eq. ^ 
disappear: 



A/(W(t)}) = J2 



Mi.jnj{t) _ iVt, 
Ntotit) 



No 



(5) 



The total population size is limited by the last term, 
which includes the carrying capacity No- The off- 
diagonal elements of the interaction matrix Af/j are ran- 
domly drawn from a uniform distribution over [—1, +1], 
while the diagonal elements are set to zero. For Model A, 
F = A and A^o = 2000 arc used in this paper. The value 
F ~ 4 for Model A is chosen such that the population 
dynamics for a single species should relax monotonically 
to a stable fixed point in order to ensure that any com- 
plex dynamics are due to interspecies interactions f20j . 
As shown in [H, [sil, [12 , communities tend to evolve 
toward mutualism in Model A. 

In Model B [l^, the external resource R is introduced. 
All the species have positive values of the birth cost 6/, 



which are randomly drawn from [0, 1], and a certain pro- 
portion (0.05 is used in this paper) of species can feed 
on the resource, i.e., the resource couplings 77/ are pos- 
itive for primary producers or autotrophs, and zero for 
consumers or heterotrophs. Here an abiotic resource R is 
introduced that is renewed each generation at the same 
level (here, 2000) and does not have independent dynam- 
ics. The off-diagonal part of the interaction matrix is 
limited to be antisymmetric. Non-zero elements are as- 
signed randomly to the pairs of (M/j, Mji) with proba- 
bility c = 0.1, which is consistent with food webs in na- 
ture, such as St. Marks Seagrass, St. Martin Island, and 
Little Rock Lake [S^HJ]. The nonzero elements of the in- 
teraction matrix are randomly chosen from a triangular 
distribution on [—1,-1-1]. These parameter ranges were 
chosen to cornparc with the corresponding individual- 
based model [1^. The diagonal elements of M, which 
represent the intraspecies interactions, arc selected ran- 
domly from a uniform distribution on [— 1 , 0] for all the 
species. The environmental carrying capacity term is not 
included in this model {No = 00). Thus, 



Ai{R, {nj{t)}) = -hi 



Ntotit) 



sr^ Mijnj{t) 



The birth cost term and the negative diagonal elements 
Mil prevent species populations from growing to infin- 
ity. The fecundity F for Model B is set to 2. With this 
value, the dynamics for a single species approaches its 
fixed point monotonically. 

These models have fixed points which can be cal- 
culated exactly [l^. Here |n*) is a column vector of the 
equilibrium population sizes of all species present in the 
community. Linear stability around this fixed point can 
also be estimated. See Appendix IB] for these solutions. 



B. Introduction of new species 

Communities are assembled through mutations of res- 
ident species as follows. Each species has a bit-string 
genome of length L, thus the total number of potential 
species is 2^. All the species-specific values, &/, rji and 
Mij, are predetermined at the beginning of the simula- 
tion and fixed during the simulation. In every genera- 
tion, a mutation happens to the genomes of the exist- 
ing offspring at the moment: all the bits existing in the 
system, which amount to iV^'ot^ x L, flip independently 
with a probability /i/L, resulting in the appearance of 
new species. The genomic mutation rate, /i, determines 
how frequently individuals mutate. Here, the number of 
individuals belonging to species /, n'/"^ is calculated by 
rounding off the population size, n'/"^ = [n/ -I- 0.5J , and 
the total number of individuals is A^tot* = J2i nf^. Thus, 
an individual moves to a neighbor in the L-dimensional 
hyper-cubic genome space by a mutation. The proba- 
bility of m-bit mutations in a single individual is small, 
0(/i™), therefore the probability of multi-gene mutation 
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is small. The coefficients of species / {bj, rji, and Mjj) 
have no correlation with those of its neighbor species, 
which is a less realistic aspect of the model. However, 
the model captures the aspect that the number of mu- 
tant species accessible from a given community is limited 
[ssl ]. Models to overcome this problem have also been 
proposed, and it is confirmed that the phenotypic corre- 
lation does not qualitatively alter the long-term fluctua- 
tions [U. 

Each simulation run was started with a single, ran- 
domly chosen species (producer species for Model B) with 
a population size of 100 individuals. The details of this 
initial condition are totally insignificant, and the systems 
were completely "thermalized" during the initial warm- 
up periods. 

We also note that the results shown below do not show 
qualitative dependence on the precise values of L and fi 
for reasonable ranges. If the mutation rate is too high, 
the system shows mutational meltdown and the number 
of species diverges. For too short L (typically L < 10), 
the system is trapped in a certain state and the species 
composition never changes. We choose parameters so 
that these unrealistic cases are excluded and simulations 
are computationally feasible. 

III. RESULTS 

In this section, we focus on the effects of the demo- 
graphic noise on the community structure. 

A. Robustness against the noise 

1. Model A 

First we show how the population fluctuations affect 
the growth of the community diversity. Figure [T] shows 
the time evolution of the diversity index D and the total 
population size A^tot for Model A with k ~ 1 and 0, 
respectively. Here, the diversity index is defined as Z) = 
exp (5), where 

Siinjit)})^- J2 Pi(t)hipiit) (7) 

{l\pi{t}>0} 

with pi{t) = ni{t)/Ntot{t), is the information-theoretical 
entropy of the species distribution. The quantity D is 
known as the exponential Shannon- Wiener diversity in- 
dex [36| . This measure can be interpreted as an effective 
number of species, and so it has the same units (species) 
as the species richness. We adopt this measure in or- 
der to filter out unsuccessful mutants which have tiny 
populations and rapidly go extinct. It has been con- 
firmed that D is approximately proportional to the num- 
ber of species constituting "core communities," composed 
of species who succeeded to have a positive stable fixed- 
population [l^. Here we do not use a finite-size cor- 
rection [s^ for the estimation of S for simplicity. This 



correction is numerically confirmed to be less than one 
percent and therefore we here adopt the simpler form 
[Eq. ©]. 

The mutation rate ^ = 0.001 and the genome length 
i = 13 are used for Model A. The difference of the aver- 
age D and iVtot between k = and 1 is slight. Both 
figures show similar intermittent behaviors, consisting 
of active and quiet periods. During the active periods, 
the diversity measure and the total population size show 
larger fluctuations, and the species composition changes 
quickly. On the other hand, during the quiet periods, 
the species composition remains nearly constant, and the 
system is considered to be in a quasi-steady state (QSS). 
The evolution proceeds intermittently, rather than grad- 
ually, like a stick-slip motion, repeating active and quiet 
periods. Average diversity and total population size for 
Model A are summarized in Table HI Both measures de- 
crease moderately with increasing noise level. 

TABLE I: Numerical results for Model A. The data are aver- 
aged over twelve independent runs. The initial 2^* generations 
are considered as a "warm-up" period and not included in the 
statistics. The statistical errors are shown in parentheses. In 
this paper, the statistical error of the value x is calculated as 
V (X] ~ i^))^) I in{n — 1)), where n is the number of inde- 
pendent runs. 



K 


D (species) 


A'tot (individuals) 


1 


3.55 (3) 


3155 (11) 


0.1 


4.81 (11) 


3246 (23) 





4.55 (9) 


3337 (36) 



2. Model B 

On the other hand, for Model B, the dependence on k 
is remarkable. Figure [2] shows typical time series of the 
diversity index and the total population size for Model B 
at several noise levels with /i = 0.0005 and L = 18. The 
data are plotted every 8192 generations for improved vis- 
ibility. As the noise level k increases, both the diversity 
and the total population size decrease remarkably. For 
other /i and L, strong dependence on the noise is also ob- 
served for Model B. Thus, Model A and Model B show 
fundamental differences. This means that it gets more 
difficult for species to survive, due to the stochastic popu- 
lation fluctuations. Although the species populations ba- 
sically fluctuate around their fixed points, the probability 
that a population size touches the extinction threshold 
increases under strong stochastic population fiuctuations. 
When K is small, high diversity and high total popula- 
tion size are realized. Especially, when k = 0, the system 
does not reach a statistically stationary state, even after 
80 million generations. In addition, the fiuctuations are 
less intermittent than for k = 1. This intermittency will 
be quantitatively estimated in the next subsection. We 
also note that the diversity for k = 1 is approximately 
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FIG. 1: (Color online) Typical time series of the exponen- 
tial Shannon- Wiener diversity and the total population size 
plotted every 8000 generations for Model A with /i = 0.001, 
L = 13, and (a) k = 1 and (b) k, — 0. The upper and lower 
curves in the figures show total population size and diversity, 
respectively. 



FIG. 2: (Color online) Time series of (a) exponential 
Shannon- Wiener diversity and (b) total population size for 
Model B with ^ — 0.0005 at several noise levels. The data 
are plotted every 8192 generations. In either figure, the curves 
correspond to k = 0, 0.1, 0.5, and 1.0 from top to bottom, 
respectively. 



the same as that of the individual-based model [1^, in- 
dicating that the current model with k = 1 is a good 
approximation to the individual-based one. The aver- 
ages of the diversity and total population size for several 
values of k are summarized in Table HH 



TABLE II: Numerical resuhs for Model B. 



The data are av- 

,24 



are considered as a "warm-up" period and not included in 
the statistics. The statistical errors are shown in parentheses. 
Uncertainties of A^tot are also in units of 100 individuals. 





D (species) 


TVtot (100 individuals) 


1 


12.2 (6) 


127 (8) 


0.5 


24.7 (4) 


164 (4) 


0.1 


101 (2) 


278 (5) 







315 (6) 



To estimate the relevance of the noise effects for Model 
B, the species abundance distributions (SADs) were in- 



vestigated. The SAD is the distribution of the popula- 
tions of each species with binning on log2 scale, which is 
widely used in ecology. Figure 3(a) shows the SAD for 
Model B with k = and 1. Since the system without de- 
mographic stochasticity does not reach a stationary state, 
we divided the time series for k = into three regions; 
ri, r2, and are the regions where t < 16, 16 < t < 64, 
and t > 64 million generations, respectively, in order to 
see how the distribution changes during the evolutionary 
process. For k = 1. the data are calculated for i > 16 
million generations, where a statistically stationary state 
is realized. During the simulation, there arc many species 
with small populations, which correspond to unsuccess- 
ful mutants. We filtered out unsuccessful mutants and 
obtained the "core fixed-point communities" by updat- 
ing the population dynamics without the noise and mu- 
tations until all the fixed-point populations analytically 
calculated from the interaction matrix (Eq. (jBip ) be- 
came larger than the extinction threshold. The SADs 
were calculated for these fixed-point communities. 
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The SADs for these communities are quite similar, 
while the peak position shows dependence on k. The 
SAD for K = 1 has a peak at higher population size than 
for K = 0. The number of species with small popula- 
tions are suppressed by the noise, and a small number of 
large-population species survive, which is a natural con- 
sequence of the law of large numbers. The difference in 
height corresponds to the difference in diversity between 
K = and 1. The profile is well fitted by a function 
suggested by Pigollotti et al. [Ill 



the population fluctuations. The distribution is shown in 



(8) 



where 7 and /3 are fitting parameters. This function in- 
terpolates between the well-known Fisher's log-series and 
the log-normal distributions. The agreement of the fit- 
ting function for k = is reasonable although there is 
some difference at n ^ 1. For k = 1, the fitting is rea- 
sonable only around the peak. The data have fatter tails 
than the function of Eq. ^ . 

Since the noise term for the fixed-point community is 
of order Ky/n/2, the noise term is of the same order as 
n* if n* < f . Simply comparing with the SADs, which 
peak at > 10^; the noises are less than the correspond- 
ing population sizes for most of the species. However, 
if the system has weak linear stability, the noise may 
have relevant effects and cause extinctions of species. The 
eigenvalues of the linear stability matrix. As, are there- 
fore estimated, and their probability density functions 
are shown in Fig. |3(b)[ The distributions for Model B 
show peaks slightly below 1, which means the system is 
asymptotically stable, but the linear stability is weak. 
The distributions show dependence on n, and the sys- 
tems with K ~ Q are less stable than those with k = 1. 
Thus, more stable communities are selected under the 
stochastic noise. 

We estimated how large the population fiuctuations 
would be if the noise corresponding to k = 1 were ap- 
plied to the fixed-point communities. Following the dis- 
cussion in [22I, we assume the probability that the sys- 
tem is found with a specific number of individuals |n) at 
a stationary state, V*{\n)), takes the Gaussian form: 



V*{\n)) = (27r)--^/2(j^^Q)-i/2 



cxp 



1 



(An|G^i|A77,) 



(9) 

where Af is the number of resident species, G is the co- 
variance matrix to be estimated, and |An) ~ \n) — |??*) 
is the difference from the fixed point. This approxima- 
tion is valid only when the population fluctuations are 
small enough to neglect the nonlincarity of the popula- 
tion dynamics. Although this assumption is not satisfied 
for Model B, this discussion tells us that the population 
fluctuations could be of the same order as the popula- 
tion sizes and might cause extinctions. How the covari- 
ance matrix G is calculated is shown in Appendix [C] 
Here we only show the distribution of the square root 
of the eigenvalues of G, which correspond to the size of 



Fig. 3(c) in the same format as the SADs (binning in log2 
scale). Comparing this figure with the SADs, the popula- 
tion fiuctuations for k = would be of the same order as 
the population sizes. Thus, the noise drives the species 
with little stability to extinction, and, as a result, cause 
the large decline in diversity. The same linear stability 
analysis was also done for Model A (shown in Appendix 
ID|) . Model A shows stronger linear stability and smaller 
eigenvalues of the covariance matrix G. Hence Model A 
is robust against the noise and does not show notable 
dependence on the stochastic noise. 

The next question is how the communities obtained 
for Model B become sensitive to the noise. The distribu- 
tion of the birth cost bj, the resource-coupling coefficient 
T]i, and the interspecies interaction coefficient Mjj are 
shown in Fig. 21 The distributions of bj and Mjj show 
dependence on k, while those of rji do not show notable 
K-dcpendcncc. In the communities which have evolved 
via population dynamics without demographic stochas- 
ticity, species with quite low bj are selected, while the 
selection on Mjj is weak. The ratio of low bj species 
increases as the evolution proceeds. On the other hand, 
under the noise, the distribution of 67 is not as extreme 
as for K ~ 0, but the ratio of species with large Mjj 
becomes larger. Hence, the selection pressure is applied 
on the birth cost for k = while it is applied on the 
interspecies couplings for k = 1. 

We also modified Model B so that all the species have 
the same birth cost, bj (= 0.1). The results are shown in 
Appcndix[El For this modified model, qualitatively simi- 
lar results as for the original Model B are obtained: large 
decline in diversity, more strongly coupled communities, 
and stronger linear stability around the fixed points are 
observed under the noise. 



B. Long-term fluctuations 

Not only the mean value of the diversity, but also its 
fluctuations during the evolution are affected by the de- 
mographic population fluctuations. The power spectral 
densities (PSDs) of the time series of diversity and total 
population size, as well as probability densities of species 
lifetimes and QSS durations, were calculated in order to 
evaluate the intermittency quantitatively. These results 
are of particular interest in connection with the dynamics 
of mass extinctions on geological time scales jsoj. 

Power laws are estimated by fitting to log2-binned den- 
sities. However, exponents obtained using the estimators 
from [40| arc not qualitatively different. 



1. Model A 

We performed simulations of 2^^ = 33554432 gener- 
ations with 2^^ ~ 4194304 generations as a "warm-up" 
period. This warm-up period is long enough to realize 
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statistically stationary states. For each model and pa- 
rameter, six independent runs were performed. Figure [5] 
shows the PSDs for several noise levels k. 

For Model A, both diversity and total population size 
generally show approximate 1// fluctuations for all val- 
ues of the noise level k. The PSDs for weak population 
fluctuations show approximate 1// power-law behavior 
over more than five decades. Thus the l/ f fluctuations 
found in the individual-based Model A [201 robustly 
reproduced, even with deterministic population updates. 
Under very strong population fluctuations, the possibil- 
ity of extinctions caused by the population fluctuations 
is not negligible, and few communities are able to per- 
sist over very many generations. As a result, the PSDs 
for high K arc not l/f like at very low frequencies. This 
effect of population fluctuations is also observed in the 
QSS duration distributio ns and the species lifetime dis- 
tributions as seen in Figs. [6(b)] and El 

Figure 6(a) shows the logarithmic derivative of the 
time series of the diversity (i.e., dS/dt), which is aver- 
aged over 16 generations. Each curve has a sharp peak 
around the center and relatively wide wings in both tails. 
The sharp peak around zero represents that the commu- 
nity is in a quiet period. The small diversity fluctua- 
tions arise from the population fluctuations of coexist- 
ing species and the repetitive emergence and extinction 
of unsuccessful mutants. On the other hand, the large 
wings represent large rearrangements of the species com- 
position. By measuring dS/dt, we can judge whether 
the system is in a quiet period where the species com- 
position remains approximately constant, or in an active 
period where the dominant species arc replaced rapidly. 
This type of profile is similar to that of the corresponding 
individual-based Model A [il,!!^]. 

The duration distributions for QSSs are shown in 
Fig. 6(b) The QSSs are estimated as the periods between 
times when \dS/dt\ exceeds a cutoff (here, 0.02). The 
distributions show approximate power laws, repro- 
ducing the result for the corresponding individual-based 
model. 

Figure[7]shows the species- lifetime distributions at sev- 
eral noise levels. The distribution is fitted by a power law 
over more than six decades, and the exponent of the ob- 
served distribution is about 2.2, which is in reasonable 
agreement with the individual- based Model A [2^. The 
distribution does not show important dependence on the 
population fluctuations. 



2. Model B 

Next we show the results for Model B. We performed 
six independent runs of 2^^ = 67108864 generations with 
2^* = 16777216 generations as a warm-up period. The 
mutation rate ^ = 0.0005 and the genome length L = 18 
were used. 

Although the species populations basically fluctuate 
around their fixed-point values, the probability that a 



species population touches the extinction threshold in- 
creases under strong stochastic population fluctuations. 
Therefore, both diversity and total population size tend 
to decrease as k increases. We also note that the diver- 
sity for K = 1 is approximately the same as that of the 
individual-based Model B [l^, indicating that the cur- 
rent model with k = 1 is a good approximation to the 
individual-based one. The fluctuations of the diversity 
for smaller k are also larger than for larger k. These 
fluctuations for small k come from introductions of mu- 
tants and extinctions of species. When k is small (e.g. 
K = 0.1 or 0), the system shows high diversity and large 
population size, which are still growing even after 80 mil- 
lion generations. In addition, the fluctuations are less 
intermittent than for k = 1. This intermittency will be 
quantitatively estimated below. 

Figure [8] shows PSDs of the diversities and the total 
population sizes at several noise levels. The PSDs of both 
diversity and total population size show power laws. The 
exponent depends on k. When k is large, the PSDs show 
approximate l/f" behavior with a = 1.3 ~ 1.4. This ex- 
ponent a is in reasonable agreement with the individual- 
based Model B [l^ although it is slightly larger. As k 
decreases, the exponent gets closer to 2; indicating that 
the diversity and the total population size both fluctu- 
ate like random walks. Thus the population fluctuations 
change not only the average value of the diversity but 
also its fluctuations on evolutionary time scales. In the 
individual-based model, the evolution proceeds intermit- 
tently, repeating quiet periods punctuated by brief ac- 
tive periods. However, in the deterministic model, the 
evolution proceeds rather gradually, and the community 
composition changes continuously. 

Since the stationary state is not realized on the time 
scale we observed for k = 0, we also calculated the PSDs 
of a corrected time series. First we calculated least- 
squares flts for the time series, and then calculated PSDs 
of the difference of the time series from the linear fit. 
The result (not shown) does not show notable differences 
from Fig. [S] 

Figure 9(a)| shows the probability density function of 
the logarithmic diversity derivative, dS/dt. Since the av- 
eraged diversity for smaller k is much higher than for 
larger k, the distribution for smaller k is quite sharp. 
This is due to the high diversity realized for weak noises. 
To eliminate this effect, the derivatives normalized by 
the average diversities arc shown in Fig. |9(b)[ We calcu- 
lated [S{t+ 16) — S{t)]/1Q X D for every 16 generations, 
where D is the average diversity. Therefore the x-axis 
of Fig. |9(b)| has the dimension of (diversity) • (time)""'^. 
The sharpness of this normalized data are almost simi- 
lar, therefore the absolute diversity fluctuations are only 
weakly K-depcndent. 

The distribution for k = 1 has a Gaussian center and 
large win gs and looks similar to the individual-based 
Model B [13. However, the distributions for smaller k 
look different. When k = 0.1, the distribution is quit e 
well fltted by a Gaussian distribution without wings [23 . 
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The distributions of QSS durations arc calculated and 
shown in Fig. [TU]for k = 1 and 0.5. The QSSs are esti- 
mated in the same way as the previous model, but dif- 
ferent thresholds are used for each k because the profiles 
show large dependence on k. The threshold for estimat- 
ing QSS is 0.1/D, where D is the average diversity. The 
distributions show approximate 1/t power laws regard- 
less of the noise level. This is consistent with the original 
individual-based Model B [l^ in which a t~^_po'wer law 
is observed in the QSS duration distribution [l^i [3 ■ For 
K = 0.1 and 0, it is impossible to estimate QSS durations 
since the distribution of logarithmic derivatives of the di- 
versity for these parameters does not have large wings. 
If we estimate the QSS with a threshold which corre- 
sponds to the Gaussian region, clear exponential decay 
is observed. Hence the small fluctuations occur randomly 
and do not have remarkable long-time correlations. 

Species-lifetime distributions for several k are shown in 
Fig. [TTJ The distributions show a reasonable flt to a 
power law for every k. The average species lifetime for 
larger k is slightly less than for smaller k, but the depen- 
dence on the noise strength is slight. Hence the lifetime 
distribution is not notably affected at the species level 
even for small k, while it is affected at the community 
level for k < 0.1. 



IV. SUMMARY AND DISCUSSION 

The effects of demographic stochasticity are explored 
for two types of biological macro-evolution models. The 
demographic stochasticity is modeled by the noise term 
of the population dynamics. For the mutualistic com- 
munities obtained by Model A, the noise does not have 
an important effect, and the deterministic description 
does not alter the picture obtained for the correspond- 
ing individual- based model. On the other hand, the 
predator-prey model (Model B) shows a remarkable de- 
cline in diversity at higher noise levels. This is because 
the deterministic population dynamics allow species to 
coexist with low linear stability, which arc pushed into 
extinction in the stochastic population dynamics. With- 
out the noise, the distribution of the birth cost bi has 
a sharp peak close to zero, while the distribution of the 
coupling constants Mjj is almost the same as the original 
distribution. With moderate noise, the selection pressure 
for small 6/ becomes less extreme, and species with larger 
Mij are selected. Hence, strongly coupled communities 
are selected under the noise, while species with low birth 
costs are selected at low noise levels. 

For model B, species must have a strong coupling to 
the external resource or strongly predatory interaction 
coefficients compared to the birth cost to sustain their 
populations. (See Eq. ([5]).) Our result that communi- 
ties with larger Mjj are selected under the noise looks 
contradictory to the classical consensus on the relation 
between stability and complexity: a community tends to 
be less stable when the interactions are dense and strong 



[4l|-|46|. This apparent contradiction is due to the anti- 
symmetric correlation of the linear stability matrix. The 
off-diagonal parts of the linear stability matrix in Model 
B is closer to an antisymmetric form. The ratio 

A/,7 _ n}{Mjj-bi) 

Ajj n*j{Mjj~bj) ^ ' 

often becomes negative since the interaction matrix M 
is antisymmetric. Eigenvalues originating from antisym- 
metric off-diagonal matrix elements arc all pure imagi- 
nary, therefore the large Mjj do not destabilize the sys- 
tem significantly. Actually, the distribution of the diag- 
onal parts A// is similar to the distribution of eigenval- 
ues, which implies that contributions of the off-diagonal 
parts are not critical. Moreover, if the average predation 
rate is large compared to the birth cost, species tend 
to have larger equilibrium populations. That makes the 
species have higher resistance against the demographic 
noise. Thus, large Mjj often leads to larger equilibrium 
populations without sacrificing the linear stability. We 
speculate that such selection of stronger coupling inter- 
actions occurs in a wide class of predator-prey population 
dynamics models under demographic noise. 

The dynamics on evolutionary time scales for Model 
B is also altered by the noise. When an appropriate 
amount of noise is applied, the system shows approxi- 
mate 1// fluctuations in the evolutionary dynamics. The 
time series consist of long quiet periods, during which 
the species compositions arc steady, and short active pe- 
riods, in which rearrangements of species compositions 
occur with relatively large-scale extinctions. The dura- 
tion distribution for the quiet periods is an approximate 
power law, and 1// fluctuations are found for the diver- 
sity index and the total population size. However, in 
the limit of no demographic stochasticity, this intermit- 
tent dynamics is replaced by a more gradual one, and 
the time series of the diversity index and the popula- 
tion size become Ornstcin-Uhlenbeck processes. As the 
noise increases, the 1//^ PSDs gradually change toward 
1// fluctuations. We speculate that this is due to the 
smallness of the mutant's population and the weak lin- 
ear stability. Since mutants are quite prone to go extinct 
under the demographic noise, QSS communities are more 
robust against the invasions of new species. Similar effect 
is also reported in another model |14| . 

The results shown in this paper indicate that models 
without noise may be remarkably different from mod- 
els with noise. Without demographic noise, communities 
with weak stability that would be destroyed under the 
noise can emerge. Since the noise effect can be more im- 
portant than suggested by a naive 1 / prediction, the 
relevance of demographic stochasticity is not limited to 
small-scale communities, such as isolated islands, lakes, 
and experimental situations in microbiology, but can also 
exist for larger-scale ecosystems. 
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Appendix A: Stochasticity in the number of 
offsprings 

It is straightforward to extend the model so that the 
number F of offspring per individual becomes random. 
Suppose the probability density function (pdf) of F, 
q{F), is approximated by a Gaussian distribution with 
mean fip and variance ap. The pdf of the number of 
individuals in the next generation, born to parents of 
species /, pi{n') for n' > nthr is then given by 



e = 



1 - (l|M-^|&) 
(l|M-i|l) ' 



(B4) 



respectively. The coefficients £ and 9 can be considered 
as an effective coupling to the external resource and an 
effective interaction strength, respectively. To find each 
n*j separately, we now only need to insert this solution 



for N,l, in Eq. dBl 

Linear stability around fixed points can also be esti- 
mated analytically. The elements of the linear stability 
matrix S are 

d{niit + l)) 
dnj{t) 

where 5ij is the Kronecker delta, and 



1 s n*j 
Ajj = 1 



M 



I J 



Rm + Y.K MiKn*K N* 



N* 



No 
(B6) 



pi{n') 



(Al) 



is the community matrix. The system is stable against 
perturbations when all the eigenvalues of S are less than 
unity in magnitude. 



where J\f[^,a'^]{x) is a Gaussian distribution with mean 
II and variance a^. Thus, the fluctuations in population 
dynamics are more enhanced when F fluctuates. The 
limit ap ^ corresponds to the model considered in the 
body of this paper. 



Appendix B: Calculation of the fixed point and the 
linear stability 

We briefly show this solution for the sake of com- 
pleteness and readers' convenience, although it is shown 
in detail in [l^ [l^. At the fixed point, the condition 
\P{R, {n*})) = 1/F is satisfied, where |P) is the column 
vector of the reproduction probabilities. Taking the log- 
arithm of this equation gives rise to M linear equations, 
where Af is the number of populated species. The solu- 
tion for In*) is 



Appendix C: Calculation of the covariance matrix G 

For the calculation of the covariance matrix G, we need 
not only the stability matrix S but also the noise matrix 
H, which is defined as the covariance matrix of the noise 
term of the population dynamics. Since the noise for 
each species is independent, H is a diagonal matrix and 
is written as 

Hi J = dijF^K^n}Pi{\n})){l - Pi{\n})) (CI) 



= 6ijK^n*j{F-l) 



(C2) 



where Sjj is a Kronecker delta, and the relation 
FPj{\nJ)) = 1 is used to derive the second equation. 
The relation between G, S. and H is 



G SGS^ = H, 



(C3) 



= M 



)R-\b)N:,,-\i){N:,,)yNo , (bi) 



where the superscript T denotes the transpose of the ma- 
trix. Although G is not simply expressed by the known 
matrices S and H. it is written in a series as 



where M^^, [77), \b), |1) are the inverse of the submatrix 
of M corresponding to the present species, and the col- 
umn vectors of 77/, bj — In (_F — 1), and ones, respectively. 
To find each Uj, we must first obtain N^^^ ( = X) ) 
follows: 



G = H + SHS^ + SSHS'^S'^ 



(C4) 



N* 



l-RS/Q {No ^00) 

where £ and Q are defined as 

(l|M-i|r7) 



£ 



(l|M-i|l) 



(B3) 



Hence, G is calculated by the following iterations: 

G, =H + SG,_iST, (C5) 

where Gq = H. We repeated this iteration until the abso- 
lute values of all the elements of the matrix (Gfc — G^-i) 
are less than 10~^. The number of iterations is typically 
of order lO''. (It depends on the species composition.) 
Since the eigenvalues of S are close to unity, many it- 
erations are necessary to obtain an accurate estimate of 
G. 
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Appendix D: Linear stability analysis for Model A 

SAD and the hnear stabihty matrix S for Model A 
were calculated. Figure 12(a) shows the SAD for Model 



A at several noise levels. The data are sampled every 
one million generations after an initial warm-up period 
of four million generations. In the same way as for Model 
B, we first removed unsuccessful mutants and obtained 
the core fixed-point communities. The profiles show little 
dependence on the noise level. We tried fitting the SADs 
by Eq. ([8]), but the fitting does not look very reasonable. 
The eigenvalue distribution of the linear stability matrix 
S and the distribution of square roots of the eigenvalues 
of the covariancc matrix G are shown in Fig. |12(b)| and 



12(b) shows that the lin- 



Fig. 12(c)[ respectively. Figure 
ear stability for Model A is much stronger than for Model 
B. Negative eigenvalues of S correspond to oscillating 
modes. The amplitude of fluctuations that would appear 



under the noise were estimated and shown in Fig. 12(c) 



The typical size of the amplitude is smaller than the typ- 
ical population sizes and, as a consequence, the resident 
species seldom suffer from the noise effects. 



Appendix E: Trial for Model B with constant birth 

cost 

The distribution of the birth cost bj for Model B at 
K = has a quite sharp peak close to zero. This emerges 
as a result of the selection and causes low linear stability. 
Since quite low birth cost is not very realistic, we also 
tried a model in which the birth cost is fixed to be 0.1 



for all species. Typical time series of diversity and total 
population size for ^ = 0.001 and L = 20 are shown in 
Fig. [T31 Averages of the diversity index D for k = 
and 1 are 219 and 19.7, respectively. Thus, there are 
significant decreases in diversity and total population size 
with increasing noise. 

Properties of long-term fluctuations are also estimated. 
Six independent simulations of 2^^ generations with 2^^ 
warm-up generations were performed. Figure [T3] shows 
the PSDs of diversity and total population sizes. For 
K = and 1, approximate 1//^ and 1/f fluctuations 
are observed respectively. This dependence on the noise 
level is similar to the original Model B. In the same way 
as the original models, the distribution of the logarith- 
mic derivative of the diversity and the duration distribu- 
tion of quiet periods are calculated. Under the noise, the 
distribution of dS/dt has a Gaussian center and wider 
wings [Fig. 15(a)| . The duration of quiet periods dis- 
tributes broadly although it is not a power-law distri- 
bution [Fig. 15(b)| . This is characteristic of the mod- 
ified model. The species-lifetime distribution shows ap- 
proximate 1 /t^ distributions regardless of the noise levels 
[Fig. [T5(c)] . The SADs [F ig. [16(a)] , the ei genvalue dis- 
tributions of S and G [Fig. 16(b)| and 16(c)| , and the dis- 
tribution of A//J [Fig. 17(b)| show the same dependence 
on K as the original Model B. The distribution of rjj has 
a different dependence on the noise level than the origi- 
nal Model B [Fig. |l7(a)| . When the noise is not applied, 
much stronger evolution pressure is applied to rji. Al- 
though there are several deviations from the results for 
the original Model B, qualitatively the same behaviors 
are observed. 
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FIG. 3: (Color online) (a) Species abundance distribution 
(SAD), (b) probability density functions (pdf) of the eigen- 
values of the linear stability matrix, As, (c) frequency of the 
eigenvalues of the covariance matrix, Ac, for Model B with 
K = 1 and ft = 0. The data are normalized in the same way 
as the SADs. For k = 0, the data are obtained for three time 
intervals. The SADs show the number of species whose pop- 
ulations are within each bin region, for each community. The 
data were sampled every one million generations, and aver- 
aged over 18 independent runs. The fitting functions Eq. (jS} 
are shown in (a) as guides to the eye. The fitting parameters 
are /3 = 2.14 and 7 = 0.0245 for k = 0, and /3 = 3.4 and 
7 = 0.0054 for K= \. Analogous data for Model A are shown 
in Fig. [12] (Appendix [D|. 
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FIG. 4: (Color online) Probability density function of (a) 
birth cost, bj, (b) coupling constant to resource, 777, and (c) 
interaction matrix elements, \Mij\, for fixed-point communi- 
ties of Model B. The data are obtained for k = 1 and k — 0. 




FIG. 5: (Color online) PSDs of (a) exponential Shannon- 
Wiener diversities and (b) total population sizes for Model A 
with /i = 0.001 at several noise levels. In both figures, lines 
corresponding to 1// are shown as guides to the eye. Data 
are averaged over six independent runs, and their statistical 
errors are also shown. 




FIG. 6: (Color online) (a) The probability density functions 
of the logarithmic derivative of the diversity, dS{t)/dt, for 
Model A with /i = 0.001 at several noise levels. The data were 
averaged over 16 generations in each run, and then averaged 
over 6 independent runs, (b) Log-log plot of the probability 
density functions of the duration of QSSs. The QSSs are esti- 
mated as the periods between times when \dS(t)/dt\ exceeds a 
cutoff (here, 0.02). The logarithmic derivative, dS(t)/dt, was 
averaged over 16 generations as in (a). The line corresponding 
to a i"^ power law is shown as a guide to the eye. 




Species lifetime (generations) 



FIG. 7: (Color online) Species-lifetime distributions plotted 
on log-log scale for Model A with fi = 0.001 and L = 13 at 
several noise levels. The line corresponding to t^^ is shown 
as a guide to the eye. 




FIG. 8: (Color online) PSDs of (a) exponential Shannon- 
Wiener diversities and (b) total population sizes for Model B 
with n = 0.0005 at several noise levels. The data are averaged 
over six independent runs and their (small) statistical errors 
are also shown. The straight lines in each figure represent 
l/f" power laws with exponents a — 1 and 2 as guides to 
the eye. The shoulder in the population-size PSD at high 
frequencies is due to self-excited population oscillations. 
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FIG. 9: (Color online) (a) Probability density functions of 
the logarithmic derivative of the diversity [S{t + 16) — S{t)] /16 
for Model B with /j. — 0.0005 at several noise levels, (b) Prob- 
ability density functions of {[S{t + 16) — S{t)]/16} x D, where 
D is the average diversity. The data are averaged over six 
independent runs and their statistical errors are also shown. 




FIG. 10: (Color online) Probability density functions of QSS 
duration for Model B with several values of k. The QSSs are 
estimated as the periods between times when \dS{t)/dt\ x D 
exceeds a cutoff (here, 0.1). The line corresponding to 1/t is 
shown as a guide to the eye. 




FIG. 11: (Color online) Species-lifetime distributions plotted 
on log-log scale for Model B with /i — 0.0005 at several noise 
levels. The line corresponding to t^^ is shown as a guide to 
the eye. 
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FIG. 12: (Color online) (a) Species abundance distribution 
(SAD), (b) probability distribution functions (pdf) of the 
eigenvalues of the linear stability matrix, As, (c) pdf of the 
eigenvalues of the covariance matrix, Ac, for Model A with 
several k. The data are sampled every one million genera- 
tions. The SADs show the number of species whose popula- 
tions are within each bin region, for each community. The 
fitting function for the SAD (k = 1) is Eq. ^ with P = 4.5 
and 7 = 0.006. Compare with Fig. E]for Model B. 
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FIG. 13: (Color online) Time series of (a) exponential 
Shannon- Wiener diversity and (b) total population size for 
the modified Model B with fi = 0.001 and L = 20 at k = 
and 1. The data are plotted every 16384 generations. Com- 
pare to Fig. [21 
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FIG. 14: (Color online) PSDs of (a) exponential Shannon- 
Wiener diversities and (b) total population sizes for the mod- 
ified Model B. In both figures, lines corresponding to 1// 
and 1//^ are shown as guides to the eye. Data are averaged 
over six independent runs, and their statistical errors are also 
shown. Compare to Fig. [S] 




FIG. 15: (Color online) (a) Probability density function of 
[S{t + 16) — S{t)]/16 X D, where D is the average diversity, for 
the modified Model B with k = and 1. A Gaussian function 
is also shown as a guide to the eye. Compare to Fig. |9(b)| 
(b) The QSS duration distribution for n — 1. The cutoff to 
detect QSS is 0.012. Compare to Fig. [TU] (c) Species lifetime 
distribution for k = and 1. A power law is also shown 
as a guide to the eye. Compare to Fig. \TT\ 
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FIG. 16: (Color online) (a) Species abundance distribution 
for the modified Model B with k — and 1. (b) The distri- 
bution of the eigenvalues of the linear stability matrix, As for 
the modified Model B with k = 1 and 0. (c) The distribution 
of the square roots of the eigenvalues of the covariance ma- 
trix, Ag. The data are shown in the same way as the SADs, 
(binning in logj scale) . Compare to Fig. [S] 
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FIG. 17: (Color online) (a) The distribution of 77/. Species 
with large rji axe more favored in the limit without demo- 
graphic noise, (b) The distribution of the off-diagonal ele- 
ments of the interaction matrix, |M/j|. Compare to Fig. 



